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Abstract. Many existing studies on pattern formation in the reaction-diffusion systems 
rely on deterministic models. However, environmental noise is often a major factor which 
leads to significant changes in the spatiotemporal dynamics. In this paper, we focus on the 
spatiotemporal patterns produced by the predator-prey model with ratio-dependent functional 
response and density dependent death rate of predator. We get the reaction-diffusion equations 
incorporating the self-diffusion terms, corresponding to random movement of the individuals 
within two dimensional habitats, into the growth equations for the prey and predator population. 
In order to have to have the noise added model, small amplitude heterogeneous perturbations to 
the linear intrinsic growth rates are introduced using uncorrelated Gaussian white noise terms. 
Eor the noise added system, we then observe spatial patterns for the parameter values lying 
outside the Turing instability region. With thorough numerical simulations we characterize 
the patterns corresponding to Turing and Turing-Hopf domain and study their dependence on 
different system parameters like noise-intensity, etc. 


1. Introduction 

Self-organizing spatial pattern formation in interacting population models is an important area 
of research to understand the distribution of the species within their habitats. Mainly reaction- 
diffusion equation models are studied analytically with supportive numerical simulations to 
understand the nature of population patches and whether they are changing with respect to 
time or not. Investigations for various types of Turing and non-Turing pattern formations in 
two dimensional prey-predator models have received significant attention [H El El HIS]. Spatio- 
temporal prey-predator models exhibit various types of self-organizing spatial patterns; either 


they are stationary or non-stationary. Stationary patterns are characterized by heterogeneous 
distribution of populations over space which do not change with time once they reach the 
non-constant steady-state. On the other hand, the non-stationary patterns change with time 
continuously and never reach any non-constant steady-state. The stationary patterns produced 
by the spatio-temporal prey-predator models for two dimensional spatial domains can be 
classified as: (i) spot patterns (cold-spots and hot-spots), (ii) labyrinthine/stripe patterns, or 
(hi) the mixture of spot and stripe patterns. These three types of patterns are mostly observed 
for parameter values taken from Turing-domain or Turing-Hopf-domain, when the diffusivity of 
two species differ by a significant magnitude. Sometimes one does find labyrinthine patterns 
as a non-stationary pattern but for parameter values well inside the Turing-Hopf domain. The 
non-stationary patterns may include periodic patterns (sometimes as periodic traveling waves), 
interacting spiral patterns and chaotic patterns. The occurrences of spatio-temporal chaos 
for parameter values lying within the Turing domain or outside the Turing-domain remains 
a debatable and controversial issue [6]. There is no unified approach in determining the 
parameter values in a systematic way or in identifying the mechanism which leads to the spatio- 
temporal chaotic patterns mm- Derivation of necessary-sufficient condition(s), either in terms 
of parametric restrictions or system characterization for the onset of spatio-temporal chaotic 
patterns for prey-predator type interacting models thus remains an open problem. 

Recently, the role of environmental noise on spatio-temporal pattern formation has been 
studied by the researchers for the noise added spatio-temporal models of interacting populations 
[El El [ig [n]. In most of the cases, the authors have been used temporally correlated and 
spatially uncorrelated coloured noise terms. Although in some investigations the alteration of 
spatial patterns due to the presence of additive coloured noise terms have been reported but 
this type of formulation is not ecologically viable as they correspond to spontaneous generation 
or inhibition of populations at the places with no population or very low concentration of 
population. On contrary, the consideration of multiplicative noise terms signify the noisy 
variation in some system parameters involved with the reaction kinetics. It is evident that 
due to the environmental fluctuations, the carrying capacity, birth and death rates, intensity 
of intra- and inter-species competition rates, predation rates exhibit random fluctuations to a 
modular extent around some average value m- Perturbation of one or more rate constants by 
the noise terms leads to spatio-temporal models with multiplicative noise terms. Several noise 
induced phenomena like noise-induced phase transitions, noise enhanced stability/instability, 
noise delayed extinction, stochastic resonance are investigated for the spatio-temporal ecological 
models considered within fluctuating environments. These analytical/numerical findings provide 
important information towards how the species establish within their respective habitats or help 
us to identify the important risk factors leading towards the extinction of species. 

The main objective of this paper is to study a simple model of predator-prey with ratio- 
dependent functional response and density dependent death rate of predator, along with the self¬ 
diffusion terms corresponding to the random movement of the individuals within two dimensional 
habitats, and study the influence of small amplitude heterogeneous perturbations to the linear 
intrinsic growth rates. 

2. The spatiotemporal model 

Consider that u = u{x^y^t) and v = v{x^y^t) represent the prey and predator population 
densities at any time t and at the spatial location (x^y) G D, where D G is a square 
domain in M with boundary dD. The spatiotemporal dynamics of prey-predator interaction 
with logistic growth for prey, ratio-dependent functional response and density dependent death 
rate for predators is governed by the following system of nonlinear coupled partial differential 
equations (PDE), in terms of dimensionless variables and parameters: 
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subjected to positive initial conditions: 
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where u is the outward drawn unit normal vector on the boundary, and a, / 3 , 7 and 5 are positive 
dimensionless parameters. The ratio of diffusivity of two species is denoted as d. The diffusion 
terms in this reaction-diffusion system for the prey and predator population are introduced to 
take care of the random movements of the individuals within two dimensional habitats. The 
zero-flux boundary condition indicates that predator-prey system is self-contained within the 
two dimensional habitat with no population flux across the boundary. 


2.1. Turing Bifurcation 

The ecologically feasible and coexisting equlibrium point for the correspoding temporal model 
are spatially homogeneous steady state for the reaction diffusion model ([l]|^ and can be obtained 
by solving the following two algebric equations, with u^v > +. 
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If we denote the coexisting equilibrium point by [u^^v^)^ then u^ and v^ satisfy the equtions 
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1 — u -= 0 =- j — ov. The explicit expression for the components of conditions 
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for their feasible existence and possible number of coexisting steady-states are discussed in 
Ref. [ 7 ] extensively. Further we assume that is coexisting homogeneous steady-state which is 
locally asymptotically stable for the temporal counterpart of the system ([l]j^, i.e., the conditions 
dll + Ci22 < 0 and 011022 — ^12^21 > 0 are simultaneously satisfied, where aii,ai2,a2i and 022 
are given by 
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The basic criteria for the diffusion-driven instability or Turing instability is that the stable 
homogeneous steady-state becomes unstable due to small amplitude heterogeneous perturbation 
around the homogeneous steady-state. Hence for the occurrence of Turing instability, first we 
need a locally stable coexisting steady-state for the corresponding temporal model. The lin¬ 
ear stability analysis for spatiotemporal model ([l][^ around the homogeneous steady-state E^ 












and the conditions required for the onset of Turing-instability are well established and used by 
several researchers. Here, we just mention the three basic mathematical criteria required for 
Turing-instability of the model ([l][^ around as follows: 


dll + Ci22 

< 

0 , 

dlia22 — ^ 12^21 

> 

0 , 

dan + Ci22 

> 

2\/~d\/and22 — < 3 - 12^21 


(7) 

( 8 ) 

(9) 


and Turing-instability sets in at the critical wave number 

2 dan + «22 ^ 

> "■ 

Satisfaction of the condition © implies homogeneous steady-state becomes unstable under small 
amplitude heterogeneous perturbation for wavenumbers k > kcr and Turing pattern emerges. 
Type(s) of emerging spatial pattern(s) (as mentioned earlier in the introduction), depend upon 
the specific choice of parametric values and magnitude of d. 

The Turing-bifurcation curve is defined by the following equation 


dan + Ci22 — 2\/~d^ aiia22 — ai2a2i> 

Solving this equation for d, we can find the critical threshold for the ratio of diffusivity dcr, 
above which the Turing patterns emerge. In the next section we are going to report the varieties 
of spatial patterns generated by the system ([l]j^ for suitable choice of the parameter values 
satisfying the Turing bifurcation conditions. 

3. Self-organizing spatial patterns 

In this section, we present the numerical simulation results and varieties of spatial patterns that 
are exhibited by the spatiotemporal model ([l|j^ for parameter values within the Turing and 
Turing-Hopf domain. If we take the parameter values a = 2.0, /3 = 1, 7 = 0.6 and 5 = 0.1 
then we find the unique coexisting equilibrium point £’*(0.2287, 0.1436) which is a homogeneous 
steady-state for the system j|l][^. £* is locally asymptotically stable for the temporal counterpart 
as the conditions 0 and Q are satisfied. In order to understand the role of diffusivity of the 
prey and the predator towards the destabilization of the homogeneous steady-state we need to 
find the quantity dcr as mentioned in the previous section. For the chosen set of parameter 
values, we find dcr — 4.9157 and hence the heterogeneous perturbations lead to Turing patterns 
for d > dcr > 

The existence and non-existence of Turing patterns are solely dependent upon the magnitude 
of the parameters involved with the model under consideration. Here we consider a and d as 
the bifurcation parameters to construct the Turing bifurcation diagram. We have presented 
the Turing bifurcation diagram in ad-parametric plane (see Fig. [^. The bifurcation diagram 
is prepared keeping /3 = 1 , 7 = 0.6 and d = 0.1 fixed, a and d are considered as the 
controlling parameters to obtain different spatiotemporal patterns. We have presented three 
bifurcation curves in the bifurcation diagram, namely Turing-bifurcation curve (blue curve), 
temporal Hopf-bifurcation curve (red dashed line) and temporal homoclinic bifurcation curve 
(black dotted curve). The coexisting equilibrium point £* for the temporal model, corresponding 
to the system ([l|j^, is stable for a < ah and it looses stability through the Hopf-bifurcation 
at ah — 2.01 (approx). The Turing instability region is the region lying above the Turing 
bifurcation curve and which is divided into two parts by the vertical line a = ah > Part of 






Figure 1. Five different types of spatial patterns observed for the parameter values lying within the Turing 
instability region (above the blue curve which is Turing bifurcation curve) and below Turing instability region. 
Choices of parameter values are marked with four different symbols based upon the resulting patterns: o cold- 
spot, □ mixture of spot-stripe, o labyrinthine, * chaotic and + interacting spiral. Hopf-bifurcation curve (shown 
by red dashed line) and temporal homoclinic bifurcation curve (shown by black dotted curve). 


the Turing domain lying in the region a > is the Turing-Hopf domain where temporal and 
spatiotemporal perturbations are both unstable. 

For the chosen set of parameter values we find only one feasible coexisting equilibrium point 

therefore we will discuss only the spatial-pattern formation around this homogeneous steady- 
state. We have obtained various patterns generated by the model under consideration through 
exhaustive numerical simulations. All numerical simulations are carried out over a 200 x 200 
lattice with time-step At = 0.01 and spatial steps Ax = Ay = 2. All numerical simulations 
are performed using the forward Euler method for the temporal part and five point finite 
difference scheme for the diffusion part. A small amplitude heterogeneous perturbations to the 
homogeneous steady-state are introduced by uq (x^, yj) = + e^ij^ vq yj) = + eryj where 

e = 0.001, ^ij and ryj {i^j = 1,2) are spatially uncorrelated Gaussian white noise terms. We 
have observed four different types of patterns within the Turing domain, cold-spot, labyrinthine, 
mixture of spot and stripe and spatiotemporal chaotic patterns. 

We have observed only cold-spot and mixture of spot-stripe patterns for parameter values 
within the Turing domain, whereas mixture of spot-stripe, labyrinthine and chaotic patterns 
are observed for parameters lying in the Turing-Hopf domain. Numerical simulations have been 
performed for systematic choices of a and d within the Turing domain in order to understand 
the entire variety of spatiotemporal patterns for parameter values within the Turing domain. 
Obtained patterns for different combinations of a and d are marked with four different symbols. 
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Figure 2. Spatial distribution of prey (i^) obtained at t = 1000 for /3 = 1 , 7 = 0 . 6 , (5 = 0 . 1 , d = 10.0 and 
(a) a = 1.95; (b) a = 2.05; (c) a = 2.15; (d) a = 2.30; and (e) d = 1.0, a = 2.15. The color bar in each panel 
represents the magnitude of u. 


The whole regime of patterns in the Turing domain is shown in the Fig. 

A sample of each type of patterns observed within the Turing domain are shown in Fig. At 
the lower panel we have presented the interacting spiral pattern which is observed for parameter 
values just outside the Turing-Hopf domain and for d — 1.0. As the Turing instability can not 
occur for d— 1.0, we can label this pattern as non-Turing pattern. To understand the stationary 
nature of the observed patterns, we have calculated the time evolution of spatial average for the 
















prey and predator densities against time corresponding to the patterns reported in Fig. and 
the results are presented in Fig. The time evolution of {u) and {v) are considered over 
a longer period of time to ensure that the final observations are far away from the transient 
behavior. Looking at the time evolution of spatial averages presented at Fig. one can 
easily understand that cold-spot, mixture of spot-stripe and labyrinthine patterns are stationary 
patterns whereas the interacting spiral pattern is a non-stationary pattern. The time evolution 
of the spatial averages are continuously oscillating and never settle down to any stationary value. 
Corresponding to the spatiotemporal chaotic pattern presented in Fig. |^d), the time evolution 
of spatial averages exhibit irregular oscillations for t > 650 and the irregularity sustains at all 
future time. Here we have presented the results up to t = 1000 from the sake of brevity. 


4. Noise Added Model 

In this section we are interested to study the effect of noise on the spatiotemporal pattern 
formation for the model For this purpose we introduce uncorrelated multiplicative white 

noise terms to the growth equations of both the prey and the predator populations, accordingly 
the noise added model is given by: 
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where and are two temporally as well as spatially uncorrelated Gaussian 

white noise terms. The averages and correlations of the noise functions are given by, 

(6 {t, X, y)) = (6 it, X, y)) = 0 

{t,x,y)^j {ti,xi,yi)) = 2aj5 {t-ti)5ix - xi)5{y - yi),j = 1,2, 

where ai and (J 2 are the intensities of the environmental driving forces. The introduction of 
the noise terms can be justified from ecological point of view also. Here we are considering 
the dimensionless model and hence we can see from ([l]j^ that the intrinsic growth rate of 
the prey is ‘F and the intrinsic death rate for the predators is ‘ 7 ’. Now we assume that the 
intrinsic growth rate of the prey and the death rate of the predator population are subjected 
to environmental fluctuations and hence we can perturb them by spatiotemporal noise terms. 
As a result, if we consider the introduction of noise terms into the said rate constants, that is, 
1 ^ 1 -h (Ji^i{t^x^y) and 7^7 + a 2 ^ 2 {t^x^y) then we find the noise added model ([n])-@ 
from the original spatiotemporal model (l]| 2 ). 

In order to understand the effect of noise on the spatiotemporal pattern formation, we have 
performed numerical simulations with the changing magnitude of noise intensities. Here we 
consider the parameter set a = 2.15, /3 = I, 7 = 0.6, 6 = O.I, d = 8.5 and assume that both 
the noise intensities are same, that is ai — a 2 — (J> The numerical simulation is performed for 
a = O.OI and the resulting distribution of the prey population is shown at Fig. Qb) along with 
the pattern observed for cr = 0 (see Fig. |^a)). There is no significant difference in the reported 
patterns as the noise intensity is small. We have also calculated the spatial averages for the 
prey and predator populations up to t = 3000 for a = 0 and a = 0.0010 and the results are 
presented at Fig. Qc)-(d). It is clear from these figures that the small noise does not change 
the system characteristics when the spatiotemporal patterns are oscillatory in nature. One can 
observe similar results for the parameter values corresponding to the stationary patterns and 
when noise intensities are not very high. 

In order to support our claim, we have prepared a bifurcation diagram for the same set of 
parameter values, as mentioned above, and considering the ratio of diffusivity as the bifurcation 
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Figure 3. Time evolution of spatial averages of the prey {u) and predator {v) populations, presented in (dark) 
blue and and (light) green colors, respectively, for t E [0,1000] corresponding to the patterns presented in Fig. 




















Figure 4. Spatial distribution of prey obtained at t = 3000 for a = 2.15, /3 = 1,7 = 0 . 6 , (5 = 0.1 and d = 8.5 (a) 
cr = 0 (b) (7 = 0.01; the color bar in each panel represents the magnitude of u. Time evolution of spatial averages 
of prey {u) and predator {v) populations in (dark) blue and and (light) green colors, respectively, for (c) (7 = 0, 
and (d) a = 0 . 01 . 


parameter. The bifurcation diagrams are prepared for the spatiotemporal model ([l]|^ as well as 
for the noise added model (lii])-@ with a = 0.01, see Fig.[^ We have performed the numerical 
simulations for d = 1.0 to d = 11.0 with increments of 0.5 and spatial averages are calculated 
for the prey population up to t = 3000. After discarding the initial transients, all local maxima 
and minima of the spatial average of the prey population are plotted against each values of ‘d’. 
If we compare the bifurcation diagrams, then it is evident that the number of local maxima 
and minima increases with the introduction of small amplitude noise terms mainly within the 
parameter regime where we find non-stationary spatial patterns in the absence of noise. One 
natural question arises whether the bifurcation diagram presented at the right panel of Fig. 
changes from one simulation to other. It does not change as we have taken the spatial average 
at each time step and the spatial domain is sufficiently large to avoid such kind of variation in 
the bifurcation diagram. 

After understanding the effect of small intensity noise on the resulting patterns, it is important 
to understand what happens with the increasing magnitude of noise intensities. Of course, we 
can find some noisy patterns when the magnitude of a is significantly large and some times it 
is difficult to identify the particular type of patterns they are representing. In order to avoid 
this type of confusion and appropriate understanding of the underlying phenomena, we have 
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Figure 5. Bifurcation diagram with respect to bifurcation parameter d with other parameters remaining same 
as in the previous figure for (a) cr = 0 and (b) a = 0 . 01 . 



Figure 6. Signal to Noise ratio measure for a — 2.15, = 1 , 7 = 0 . 6 , 5 — 0 . 1 , d = 8.5. 

calculated the signal-to-noise-ratio (SNR) against a range of noise intensities. Keeping a = 2.15, 
/3 = 1, 7 = 0.6, 5 = 0.1, d = 8.5 fixed, we have performed the numerical simulations for a range 
of values of a. After discarding the initial transients, we have calculated the signal-to-noise-ratio 
measure obtained at t = 3000 and the average is taken over the entire 

spatial domain m- The plot of the SNR measure is shown at Fig. when a varies within the 
range [0.0025,0.1]. From the figure, it is evident that the relative variance of the population is 
not a monotonic increasing function of the noise intensities. The environmental driving force 
can cause the dispersion of the populations from their average value but the internal mechanism 
can compensate to maintain the size of the population patches and as a result we find that the 
relative variance is an increasing function of the noise intensity. 

5. Conclusion 

In this paper, we have studied a predator-prey model with ratio-dependent functional response 
and density dependent death rate of predator, incorporated with the self-diffusion terms 
corresponding to the random movement of the individuals within two dimensional habitats. 
Extensive numerical simulations were performed to understand the regimes of patterns that the 
model exhibited within and outside of the Turing boundary. Typically, four types of patterns 
were observed within the Turing domain and interacting spiral pattern as non-Turing pattern. 













The effect of environmental driving forces (noise) were also studied. Interestingly, we observed 
that the small amplitude noise does not alter the stationary pattern; only the time taken to 
reach the stationary pattern increases. However, irregularity increases even with small noise 
intensities within the non-stationary regime. We observed that dispersion of population from 
their spatial average does not increase gradually with the noise intensities; rather noise and 
demographic interaction play crucial roles to determine the distribution of populations within 
their habitats, even during the adverse environmental conditions. Our next attempt will be 
to address the issues of parameter values within the spatiotemporal chaotic domain, and also 
investigate different noise scenarios such as ai ^ (72. 
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